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^^ , We present the derivation of an interatomic potential for the iron phosphorus system based pri- 

(~| ■ marily on ab initio data. Transferrability in this system is extremely problematic, and the potential 

'—{ ' is intended specifically to address the problem of radiation damage and point defects in iron con- 

taining low concentrations of phosphorus atoms. Some preliminary molecular dynamics calculations 

VQ ■ show that P strongly affects point defect migration. 
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C/3 : I. INTRODUCTION 

Phosphorus is one of the major causes of embrittlement of nuclear reactor pressure vessel (RPV) steels[l|. A huge 
enhancement of the concentration of phosphorus (P) atoms at grain boundaries is observed in samples of RPV steels 
taken from nuclear reactors J3|- This leads to a decrease in the grain boundary cohesion and consequently to a shift 
C^ 1 of the ductile-to-brittle transition temperature. 

A theoretical understanding of the problem requires a full understanding of the interaction of millions of these two 
i_j , atoms at the atomic level. First principles calculations provide the most reliable way to describe interactions, but 
they are impractical for large scale molecular dynamics. Hence there is a need for accurate descriptions of the energy 
based on interatomic potentials which do not treat the electrons explicitly. 
fi • A crucial aspect of atomic- level modelling is proper relaxation of atoms around the defect, for which elegant methods 

—^ ] were developed by Michael NorgettQ • With state-of-the-art ab initio methods it is now possible to treat a few hundred 

^ • movable atoms: just as with empirical potentials in the seventies this brings problems of finite size effects to the fore 

and many of the concepts pioneered by Michael Norgett and encapsulated in his DEVIL 4] (Defect Evaluation In 
Lattices) code are being revisited today, and indeed many of the applications such as dislocation core structure |5j and 
V/-^ ' twin boundariesjg are still debated. 
C^ . Here we present the derivation of an interatomic potential for the Fe-P which can be evaluated as quickly as a short 

^^D ' ranged pairwise potential. For molecular dynamics purposes in studying reactor steels, the interesting region is that of 
^^ '\ small concentrations (~ 10"'^) of P in Fe, in particular, the behaviour of point defects in lattices. Traditionally, data 
for fitting potentials came from experiment, typically bulk properties, but recently ah initio total energy calculation 
enables us to add more specific atomistic level configurations to the fitting dataset. Here we incorporate both ab 
t^ ' initio and experimental data to parameterise a potential for the dilute P-Fe system. 
We choose to write the potential in the form: 
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For Fi[x\ = ^/x this is the second- moment tight-binding form of Cyrot'7'l and Finnis and Sinclair 8]. This is 
transferable between environments where the local band structure is primarily changed by scaling 9] . Wc will need to 
fit seven functions for the two-component system: three pair potentials VFeFeif), VFepif) and Vpp{r), three pairwise 

^ . functions (fiFeFeir), (t^Fepir) and (f>pp{r), and two embedding energy functions Fpe and Fp. 

■ ■ ' A previous parameterization in this form for general properties of iron|lClj gave good qualitative results but was 

not tailored specifically to defect properties. Recent ab initio work[l3[l3 shows that at high densities the resistance 
to compression arises from a many-body effect (electronic kinetic energy) rather than a pairwise repulsion. Fitting 
compression data using the short range part of V{rij) typically leads to an overestimate of interstitial formation 
energies and volumes 14]. These quantities have recently come into the realm of what can be calculated by ab initio 
means, and can be used to fit short range repulsion in regions where the density is close to its bulk value. We have 
shown that including different things in the fit leads to very different potentials: to describe point defect interactions, 
one must fit poi nt defect properties, and so we use an iron potential optimised for point defects as the basis of the 
present work[l3. This incorporates the kinetic energy effect by writing F{x) = —^px^a-ix^ -1-040;^, with increasing x 



taking us smoothly between hopping and free-electron dominated regimes. Hence isotropic compression is dealt with 
by the many-body term (high x), while short atom-atom distances can be addressed with the pair potential V{rij). 
This modest change to the second-moment formalism gives greater fitting flexibility and no additional computing 
cost, since F is implemented as a look up table. 

Pure P is covalently bonded and cannot be described by this type of potential. We therefore do not attempt to 
fit properties of pure phosphorus, or phosphorus-rich compounds, concentrating instead on point defects in a— iron, 
their interactions with phosphorus atoms and fictitious iron-rich compounds. To parameterise the potential we use 
data generated from ab initio plane wave pseudopotential calculations using density functional theory 16], ultrasoft 
pseudopotentials[i3, the spin-dependent generalised gradient approximation for exchange and correlation[i^ and 
periodic boundary conditions. Calculations on pure iron[l^ suggest that this is a reliable approach, it has been 
deployed widely and the calculations can be routinely done using standardised software[20J. 

There is an additional caveat for ferromagnetic materials: The magnetism (and hence Fermi level) is affected by 
defects and this leads to a much slower convergence of the energy with system size (and k-point sampling) than is 
typically observed for non-magnetic elements |23|- This may be due to the fact that the Fermi energy moves relative 
to the spin-bands in the supercell calculation, whereas for a truly infinite crystal the bands are fixed. Despite this 
slow convergence of defect energy, the energy differences between e.g. different interstitial configurations converge 
much more rapidly[l9j. Thus we are justified in fitting energy differences from calculations using relatively small unit 
cells, while obtaining the absolute formation energy from larger calculations. 

Ferromametism also means that the band structure changes dramatically at the fcc-bcc phase transition iron. For 
this reasonQ we do not expect the potential to describe paramagnetic iron correctly: high-T fcc-bcc transitions have 
been observed with other iron potentials but the driving force is vibrational (entropic) not magnetic |15J |. 

II. METHODS 

The fitting strategy assumes that phosphorus interacts with point defects via short-ranged, pairwise interactions 
and long-range strain fields. Thus we include configurations representing point defects, single substitutional impurity 
(SSI) phosphorus atoms, and combinations thereof. Since the major problem associated with phosphorus is segregation 
to grain boundaries, we include Fe with a 2D layer of P in our fitting database. By fitting relaxation volumes, we 
capture long-r ang e strain effects. We also include liquid configurations to ensure that the functions are sampled at 
all separations |3l|. 

The ab initio calculations and configurations included in the fitting are given in table JlHllIIII Our ab initio 
calculations use small unit cells which introduces problems of images and relaxation:4] for considering isolated defects. 
However, by fitting the same small cells described by the potential, we alleviate this problem. 

The most crucial aspect for radiation damage is the geometry of the defects and barriers (which governs diffusion 
mechanisms) and the energy differences between them (which governs diffusion rates). "Geometry" here includes the 
symmetry of the defects and their P-composition, but not detailed interatomic separations. Similarly important are 
the interaction strengths between phosphorus and point defects, which determine whether pinning of defects occurs. 
Of secondary importance are the actual values of the formation energies: this affects production rates in cascade 
simulations, but in a molecular dynamics run defects are very unlikely to be generated thermally. 

A. Ab initio Calculations 

Our ah initio calculations are done using standard codes |20j, implementing the pseudopotential plane wave method 
using density functional theory'lG^, ultrasoft pseudopotentials'lT^ , the spin-dependent generalised gradient approxi- 
mation for exchange and correlation 18] and periodic boundary conditions. The plane wave cutoff of 300eV has been 
used throughout with fc-point convergence to O.OleV, force convergence to O.OOleV/A and strees to O.OlkB. 

Previous ab initio calculations on iron-phosphorus |22. 123 . |24L l25l ] using the similarly-reliable full-potential linearized 
augmented plane wave method have downplayed total energy and concentrated on the role of increased electron density 
as indicating strengthening of the bonds at the surface relative to the grain boundary. This qualitative picture is not 
specific to FeP and follows from the tight-binding formalism[3. We observe similar increased electron densities in our 
calculations of a single layer of P in Fe, but do not use them in fitting. 

The actual energies calculated by DFT relative to the free atom are known to be unreliable. Likewise, the Finnis 
Sinclair formalism cannot be expected to be transferrable to free atoms. Thus we do not at any stage fit the ab 
initio energies directly, rather, we fit relative energies between different configurations of the same number and type of 
atoms. We take the experimental cohesive energy for a-iron (4.316cV/atom), but the cohesive energy for phosphorus 



is not fitted (and nor is the crystal structure of pure phosphorus). All fits are to energy differences with phosphorus 
in various locations, with the single substitutional impurity taken as the reference state. 

1. Crystal structures - (table 1) 

Several crystal structures were examined and used in the fitting. At the ¥c^V composition we looked at fictitious 
LI2 (fee equivalent), DO3 and DO32 (bee equivalents) and a bcc-based structure with every fourth (001) layer replaced 
by phosphorus. DO3 was noticably lower in energy. 

We also calculated Fe2P, a complex structure which does exist experimentally, finding excellent agreement a=5.836 
(expt 5.865), c=3.436 (expt 3.456), u=0.257 (expt 0.256), v=0.591 (expt 0.594). This gives us confidence that the ab 
initio calculations are reliable for the system. Results are given in TableQ] 

2. Monolayers and surfaces (table 1) 

To include fracture-relevant data, we evaluated surface energies for pure Fe, Fe with a monolayer of P at the 
surface, and Fe with a monolayer of P one layer below the surface (there is evidence that P-rich grain boundary 
fracture occurs by breaking of Fe-Fe bonds adjacent to the boundary, rather than Fe-P bonds at the boundary.) The 
results show that phosphorus does not segregate into a monolayer, either in bulk or on the surface, and that a free 
(001) surface with P on it has a higher energy than without (taking substitutional phosphorus as the reference state). 
Assuming replacement of iron by P, the segregation to grain boundaries must result from the different crystalline 
environment there rather than an intrinsic tendency of P to form layers (by extension, some grain boundaries will 
be more susceptible than others). Moreover, the fracture embrittlement probably arises from a stress concentration 
effect rather than a simple Griffith-criterion energy balance. We note however that some experimental evidence 2^ 
suggests that P sits in hollow sites on Fe(OOl) - investigating all such possible reconstructions by ah initio calculation 
would be impractical, hence the need for reliable potentials. 

3. Substitutional impurities (table 1) 

Phosphorus has a small range of solid solubility in Fe, and is believed to be located substitutionally. Several 
configurations were calculated: a single P atom in a 16(54) (2'^(3^) bcc) supercell and a pair of P atoms on neighbouring 
sites in similar cells. In each case all ions were relaxed and we took the calculations to k-point convergence (729/216 
k-points). We compared the effects of fixing the lattice parameter at the pure Fe value, or relaxing it to minimise 
energy (see Table. EJ. The finite size effects going from 16 to 54 atoms were only about O.leV, within the fitting 
errors. In practice the FeisP-SSI value is used as the P reference state. 

The addition of P impurities causes a striking reduction of 10-15% in the stiffness of iron. We calculate elastic 
constants by applying finite strains and measuring resultant stress. Our calculations for pure iron give elastic constants 
of Cii=225GPa, Ci2=124GPa, C44=101GPa, some 10% lower than experiment (and therefore not included in the 
fitting). Calculations on the FeisP substitutional impurity supercell (6.25% P) give Cii=196GPa, Ci2 = 109GPa, 
C44=91GPa. 

4. Vacancies (table 2) 

With a 31-atom Fe unit cell, we find a vacancy formation energy of 1.94eV compared to the fully converged value 
of 1.95eV[i3. With 30-1-1 atoms, the energy required to create a vacancy adjacent to a substitutional phosphorus is 
1.64eV for the near neighbour site and 1.72eV for the second neighbour site. 

For 15-atom supercells, the energies are 1.71eV for the near-neighbour site and 2.18eV for the second neighbour 
site, this latter suggests a P-vacancy repulsion, but the cell is so small that the calculation actually represents an 
alternating chain P-vacancy-P-vacancy, so we neglect it. 

We estimate the migration barrier by replacing two atoms in pure Fe with one located at their midpoint [23 and 
relaxing the remaining atoms while constraining the mirror plane symmetry. As we shall see, in bcc this is not 
necessarily the barrier for a (^, i, i) hop but it does represent a useful configuration to include in the fitting. The 
energy for P to hop from the near neighbour site into the vacancy is much lower than for the direct second neighbour 
[001] hop: in the fitting we ensure that the latter barrier is high enough not to be surmountable in MD. 



5. Interstitials (table 2) 

In common with experiment |23 and previous work |10lll9l | we find the [011]Fe-Fe dumbbell configuration to be the 
stable interstitial. Geometrically, the migration mechanism for an [Oil] interstitial could go via the [III], tetrahedral, 
octahedral or [001] configuration - our calculations suggest the tetrahedral and [III] configurations are similar in 
energy. Dynamically, it is also possible that a process of excitation to (III) followed by fast ID migration may be 
favoured|lO|. 

Although we use smaller supercells without relaxation, the difference in energies from larger calculations can be 
approximated by elasticity [23, and the differences between various conformations (which is what we fit) are converged 
to within the accuracy with which we can fit them. When one or two phosphorus atoms form part of the interstitial 
the energy is lowered compared with an iron interstitial and substitutional phosphorus, and the stable configuration 
remains the [Oil] dumbbell. 

6. Liquid state calculations 

We use a self-consistent process to model liquids, starting with pair potentials|23 and MD to create a "typical" 
atomic-level model of liquid Fe-P alloy (84 Fe atoms, 12 P). Then the total forces acti ng o n each atom in this model 
are obtained from static first principles calculations. We then use dynamical refitting[I2 via force matching[33 to 
produce a new trial potential and generate a new "typical" liquid configuration with classical MD. The process is 
then iterated to self-consistency. Liquid configurations provide data across the range of possible Fe-P separations, 
including small separations which are absent in equilibrium crystal data at T=0, and ensure there are no anomalous 
wiggles in the potential at separations for which there is no data. 

B. Fitting 

In addition to the ab initio calculations above, other properties of pure iron fitted to experiment are 
Cii=I.5I7eV/A3, Ci2=0.86IeV/A3, C44=0.76IeV/A3, and cohesive energy of 4.3I6eV. 

We choose to parameterise the potential using a polynomial spline functional form. This choice is arbitrary, but 
does not constrain the physical behaviour - we have shown 31] that a liquid simulation using a cubic spline fit can be 
used to fit a reciprocal power series potential, and the resulting V{r) and (pir) are indistinguishable. 

Thus 

<^(r,,) = Y. MRk - rfH{Rk ~ r) (2) 



F{x) — ~\fx + a-ix^ + a^x (3) 

For the extreme short range repulsion, which is sampled only by the primary knock-on atom in a cascade simulation 
and not by any of our fitting data, we adopt the screened electrostatic form of Biersack and Ziegler 32]. 
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+ H{r2 - r)H{r - ri) exp (Bq + Bir + Bar^ + B^r^) 

+ Hin-r)^ar/rs) (4) 

where Qi and Qj are the nuclear charges and where Vg = 0.4683766/(Qj -I- Q ). 

C(x) = O.I8I8e-3-2^ -I- 0.5099e-"-^^23^ + 0.2802e-"''"29^ + 0.028I7e-°-2°i6^ (5) 

In these equations, H{x) is the Heaviside function and rs the Bohr radius. The Bi coefficients are determined by 
continuity of the potential and its derivative, and so the parameters available for fitting are the Ak, au and C. 

We take Vp^Fe^ 4'FeFe and Fp^ from our previous paper [l^. The (jis are taken by scaling the pure iron values. 
The implicit assumption is that in the rigid band picture, the tight binding energy goes parabolically with number of 



valence d-electrons: (lO-N)N . We assume that in the dilute alloy the 3 phosphorus p-electrons end up in an unchanged 
d-band. This suggests that the bond strength is e.g. '^Fep/'^FeFe = (3(10 — 3)/6(10 — 6))^ = 0.765625 the squaring 
coming from the assumption that the embedding function is a square root. Obviously, this means that the potential 
is invalid for high phosphorus concentration. 

To fit the remaining parameters, we use a weighted least squares fit to the lattice parameter and formation energy 
for several FesP compounds (DO3, LI2, stripe), the relative energies of defect structures (including the (100), (110) 
and (111) mixed interstitial configurations), vacancy -substitutional interaction and vacancy migration energies from 
the first principles calculations. 

In the weighting, prime importance was placed on those configurations expected to be seen in simulation, the 
stable [110] interstitial, the substitutional phosphorus and the migration barriers. While other configurations were 
less strongly weighted, we ensured that they remained sufficiently high in energy that they would not subsequently 
participate spuriously in dynamics. 

For the defects, as with the liquid, the fitting process is a self-consistent one - we fit to unrelaxed defect configurations 
to obtain a trial potential, relax the configuration using this potential, then refit the potential to the new configuration. 
After several iterations, a self-consistent set of defect and liquid configurations, energies and potential parameters is 
obtained. This process means that the interatomic distances found in the ah initio are not fitted. Once a fit had 
been obtained, the functions V{r) and 4>{r) were tested ab oculo for overfitting and MD simulations for thermal 
expansion, elastic moduli (FigQ, vacancy and interstitial migrations and liquid state diffusion were done as a check 
for pathologies. The parameters for the final potential are given in table IVl 

III. PHOSPHORUS DIFFUSION MECHANISMS 

A. Interstitial Mechanism of Migration deduced from ab initio 

The models of P atom segregation proposed so far jsi'', 35, 36] are developed for the case when the binding energy of 
a phosphorus-interstitial complex is small. Specifically that the mean free path until thermally-activated dissociation 
is much shorter than the mean distance between sinks of point defects. This condition is needed to justify the 
detailed balance approximation, required for the concentration of P-interstitial complexes to be a function of the local 
concentration of P atoms. Our ab initio calculations show that the binding energy of a Fe-P mixed dumbbell is not 
small, and hence this approximation is not valid. Indeed, the large binding energy and small migration energy of 
an Fc-P interstitial complex implies that during the lifetime of this complex the probability that it will dissociate 
thermally is small. For a typical value of phosphorus concentration in RPV steels, ~ 10^'^, the mean free path of 
an irradiation-produced Fe-Fe dumbbell before encountering a P atom is several nanometres. Hence, the majority of 
interstitials arriving at a GB (or any other sink of point defects) will contain phosphorus atoms. As a consequence of 
the phosphorus-interstitial complex should be treated as a single migrating entity in methods such as kinetic Monte 
Carlo. 

B. Vacancy Mechanism of Migration 

Our results also show significant vacancy-P atom interaction energy: the binding energy of this complex to be 
~ 0.3eV. This was not expected previously and in the model by Lidiard|33 this interaction was totally neglected. 
Even more importantly the interaction is long-ranged, at least up to the second neighbour atoms, which invalidates the 
simple models 34] for diffusion coefficients in Fe-P alloy. Their conclusion, that in bcc alloys solute atoms always drift 
up the vacancy concentration gradient 34] , must be re-examined allowing for longer ranged interactions and complex 
formation. The long-range vacancy-phosphorus interaction in bcc iron makes it possible for a vacancy to move around 
the P atom, while remaining bound as a complex. As a consequence the situation may become similar to fee alloys, 
where the diffusion coefficient of a solute atom can be positive (indicating drag of solute atoms) or negative (simple 
exchange) depending on the relative frequencies of vacancy jumps between two neighbour sites of the solute atom and 
away from the solute atom. Hence, it is possible that vacancies would drag P atoms to sinks of point defects. 

The small energy for the vacancy P atom exchange jump (relative to that of vacancy iron atom) implies that the 
diffusion coefficient of phosphorus atoms via a vacancy mechanism would be independent of this energy. The potential 
shows that the barrier for first-to-second neighbour jumps is similar to barriers in pure Fe. Finally, the conclusion |35| 
that the migration of phosphorus atoms to grain boundaries is predominantly via the interstitials has to be verified 
in the view of strong long ranged interaction of vacancy P atom complexes, and smaller production rates of single 
interstitials than single vacancies in high-energy displacement cascades in a-iron due to higher intra-cascade clustering 
of interstitial atoms. 



C. Calculated Mechanism of Migration using the Potential 

We have performed simulations of diffusion for point defects in FeP using the potential defined in tabled We 
find that vacancies diffuse freely in iron, but are attracted to and form complexes with the P atoms. Although 
the phosphorus atom can move into the vacant site, in isolation this would simply produce a thermally activated 
oscillation, and no net diffusion - the vacancy has to hop out to second neighbour sites and back again for comigration 
to occur. Using static relaxation, we calculated the various barriers to vacancy jumps in the vicinity of a P atom 
(FigI21). A near-neighbour hop in bcc crosses two intermediate (111) planes, so the barriers tend to be bimodal. 
Although the P-v migration barrier is only half that of the Fe-v barrier, to obtain long range P migration the vacancy 
has to hop around the P via second neighbours. Energetics (both ab initio and with the potential) show the P to be 
bound to the vacancy in both first and second neighbour sites, however from the second neighbour site the barriers 
are similar for hopping back or further away from the P. Thus P acts as a strong vacancy trap, and comigration is 
likely. 

We performed some preliminary MD simulations using the potential with 2000 atoms to check the mechanism 
and ensure against pathological behaviour. For 1200K, with a single P interstitial, a series of ID migration steps 
were observed for about 80ps, after which the interstitial dissociated from the P atom and standard 3D migration 
in pure Fe continued. At 920K and 600K, however, the mixed interstitial did not separate and migration was more 
three dimensional (Fig|31). Interstitials diffuse freely in pure iron, and are strongly attracted to substitutional P. 
The P-Fe mixed dumbbell is also highly mobile, and the P diffuses rapidly. Thus P-diffusion occurs primarily via 
interstitials in radiation damage conditions where Fe-Fe interstitials are commonplace. Under thermal conditions 
however, interstitial mediated segregation to boundaries will not occur since there is no source for interstitials in the 
bulk, and mixed interstitials formed at sinks cannot transport further P-atoms to the sink. 

The stability of small complexes suggests, however, that still larger complexes may be important. For example 
mixed interstitials may be pinned by other P atoms, forming sessile PP interstitials. Future molecular dynamics can 
incorporate features neglected by this model, such as clustering, recombination, segregation and the specific defect 
distributions associated with thermal, electron or neutron radiation. 

IV. CONCLUSIONS 

We have developed an interatomic potential for the a-Fe-P system in the limit of low-P concentration with particular 
importance attached to configurations observed under radiation damage. This potential gives a good description of 
point defect properties and conformations, in particular vacancy and interstitial complexes with a P atom and hence 
can be used in larger scale molecular dynamics or Monte Carlo simulations to study diffusion properties and segregation 
of P during both ageing and irradiation conditions. Pure phosphorus cannot be described by this type of potential. 

The potential predicts that phosphorus is strongly bound to both vacancy and interstitial defects, in accordance 
with ab initio results. This fact has not been expected previously, and conclusions of theoretical studies that ignore 
strong P-defect binding should be reexamined. 

We have presented prelimnary molecular dynamics studies of P diffusion in iron which show the potential to be free 
from pathologies, but a quantative study of these processes (which should include larger complexes) in beyond the 
scope of this paper. It appears that in the temperature region of interest for RPV steels P would migrate rapidly via 
both interstitial and vacancy complex mechanisms, although we have not studied large complexes which may further 
complicate the process. 

In summary, the detailed understanding of radiation damage in steels remains an area of active interest more than 
thirty years after the pioneering work of Michael Norgett helped pose the questions. 
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Structure 


Energy (eV) 


Mag. Mom 


Volume A^ 


Pressure 


Fea 


bcc 


16.617 


4.62 


23.24 





FeaP 


LI2 


30.755 


7.85 


47.59 
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DO3 


30.971 


5.46 


43.55 





FegP 


stripe 


30.786 


5.52 


44.69 





FeePa 


FeaP 


69.718 


9.13 


101.34 





FcisP 


SSI 


130.867 


35.30 


185.9 


1.34 


FesaP 


SSI 


446.714 


125.26 


630.4 





Fei4P2 


DSI 


128.407 


32.71 


185.9 


6.67 


Formula 


Strain 


Energy (eV) 


Mag. Mom 


Volume A^ 


Stresses 


FctP 


ezz = -0.025 


63.997 


14.6 


88.24 


65, 82 


FerP 


ezz = 0.0 


64.066 


14.8 


90.51 


25, 10 


FerP 


ezz = 0.025 


64.042 


15.0 


92.77 


-13, -51 


FerP 


tzz = 0.050 


63.941 


15.2 


95.03 


-45, -98 


FerP 


ezz = 0.075 


63.789 


15.6 


97.30 


-70 -130 


FerP 


ezz = 0.5 


61.073 


10.9 


135.76 


-26, 


FeePFe 


ezz = 0.5 


60.988 


6.9 


135.76 


-24, 


Feg 


ezz = 0.5 


64.141 


21.0 


135.76 


-6, 



TABLE I: Energies of relaxed iron-rich crystal structures ("stripe" is a layer bcc structure with Fe-Fe-Fe-P (001) layers), single 
substitutional impurity phosphorus atoms (SSI), near-neighbour double substitutional impurity (DSI), and an (001) monolayer 
of P in 7 layers of Fe, (labelled FerP) both under strain and after fracture at the P (FeyP), one layer above the P (FegPFe) and 
in pure iron (Feg). Stresses are parallel and perpendicular to the long (z) axis respectively. The quoted energies are relative to 
a spin-zero GGA representation of an iron atom. Absolute values are not used in the fitting, only differences between them for 
which this arbitrary zero cancels out. 



Formula 


Structure 


Energy (eV) 


Mag. Mom 


Volume(A^) 


Pressure (kB) 


adjusted 


energy (eV)[22] 


Fei5P2 


111 


132.470/4.636 


33.8 


185.9 


218 


2.91 




Fei5P2 


001 


132.464/4.642 


33.3 


185.9 


246 


2.45 




Fei5P2 


Oil 


133.419/3.687 


34.1 


185.9 


223 


1.88 




FeisP 


TET 


134.828/4.347 


32.8 


185.9 


203 


2.85 




FeieP 


001 


134.428/4.747 


35.8 


185.9 


252 


2.45 




FeieP 


Oil 


135.353/3.822 


32.7 


185.9 


198 


2.39 




FeieP 


111 


134.820/4.355 


35.7 


185.9 


224 


2.54 




FeieP 


OCT 


134.428/4.747 


35.8 


185.9 


252 


2.45 




Fei7 


TET 


135.723/5.521 


31.4 


185.9 


184 


4.29 




Fei7 


Oil 


136.419/4.825 


34.2 


185.9 


201 


3.36 




Fei7 


001 


134.779/6.466 


34.6 


185.9 


224 


4.65 




Fei7 


111 


135.746/5.498 


34.9 


185.9 


201 


4.03 




Fei7 


OCT 


134.530/6.715 


35.0 


185.9 


228 


4.83 





TABLE II: Calculated results for interstitial formation energies, relative to GGA free atoms / pure iron and SSI phosphorus. 
Fei7 denotes iron interstitials, FeieP denotes mixed dumbbells or P interstitials, Fei5P2 denotes P-P dumbbells. TET and 
OCT denote tetrahedral and octahedral sites respectively, other structures are dumbbell configurations. In each case volume is 
constrained to the 16-atom pure iron cell, a 9x9x9 k-point grid was used and all atomic positions relaxed. The mixed dumbells 
have low symmetry and may relax to a higher symmetry state: [001] mixed dumbell relaxes to the octahedral position and the 
[111] mixed dumbell goes to the crowdion position. These calculations are for small supercells incorporating both interstitial 
formation energy and strain interactions - they should not be regarded as energies of isolated defects (they are probably upper 
bounds). Cohesive energies are quoted firstly relative to non- magnetic free atoms with GGA secondly relative to solid Fe and 
substitutional FeisP. Elastic correction for finite size 29] of P^V/2B suggests that full relaxation in a large supercell would 
lower the formation energies systematically by about 1.5eV 



Formula 


Structure 


Energy (eV) 


Mag. Mom 


Volume (A=^) 


Pressure(kB) 


Feis 


vacancy 


122.606/2.021 


36.4 


183.66 





Feis 


bar vacancy 


122.008/2.619 


36.4 


184.56 





Fei4 


divac-lst 


112.563/3.757 


34.9 


180.20 





Fei4 


divac-2nd 


112.838/3.481 


35.5 


179.10 





Fei4P 


1st 


120.852/1.706 


33.2 


181.78 





Fei4P 


2nd 


120.379/2.180 


34.6 


184.64 





Fei4P 


bar-Ill 


120.504/2.055 


33.7 


184.14 





FeaoP 


1st 


253.873/1.622 


71.1 


371.87 


-15 


FeaoP 


2nd 


253.775/1.719 


72.3 


371.87 


-7 


FeaoP 


bar- 100 


252.006/3.489 


70.8 


371.87 


-19 


FeaoP 


bar-Ill 


253.525/1.970 


71.1 


371.87 


-8 



TABLE III: Calculated results used to fit vacancy formation energies, relative to GGA free atoms/pure iron and SSI phosphorus. 
Formulae give number of atoms in the supercell, "divac" represents two vacancies at nearest and second neighbour sites 1st 
and 2nd represent the site of the P relative to the vacancy, "bar" represent the energy with the P midway along its migration 
path. In 15 atom cases unit cell and atomic positions are fully relaxed, and a 9x9x9 k-point grid used. These vacancy-rich 
configurations incorporate both vacancy formation energy and strain interactions - they should not be regarded as energies of 
isolated defects. In particular, for the second neighbour cases we have a chain of second neighbour vacancies. Comparison with 
results for larger supercells [iSll suggests that full relaxation in a large supercell would lower the energies slightly, but maintain 
the differences between them. 
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Structure 


Energy (eV) 


Dilatation(%) 


pure iron 


4.013 





Fe vacancy 


1.71 


-22.3 


Divac-lst 


3.285 


-35.3 


Divac-2nd 


3.18 


-48.9 


Bar- Vacancy- Fe 


2.34 


- 


Int 110 


3.59 


4-124.7 


Int OCT 


4.22 


+102.1 


Int 100 


OCT 




Int 111 


110 




PSSI 


0.0 


-34.3 


P-vac (1st) 


1.34 


-57.1 


P-vac (2nd) 


1.37 


-58.1 


Bar P-vac 


1.65 


- 


P-110 


2.57 


4-115.7 


P-111 


3.30 


+102.7 


P-100 


Int 110 




P-TET 


2.80 


+151.5 


P-OCT 


3.47 


+161.5 


Bar P-SSI-110 


0.27 


- 



TABLE IV: Energies and dilatations of defect configurations calculated using the interatomic potential, with 2000 atom constant 
(zero), pressure static relaxation. Notation is as for previous tables. Energy refers to a reference state with an equivalent 
number of iron atoms in pure iron and phosphorus atoms as SSIs in iron. For unstable interstitial configurations we indicate 
the local minimum. Dilatation is defined relative to pure iron and calculated using constant volume and elastic constantS;37l| 
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FIG. 1: Effect of 10% P on Fe lattice parameter (squares) and bulk modulus(circles), evaluated by molecular dynamics on 
2000 atom cell. For comparison similar quantities for pure Fe (red lines) are shown. 
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FIG. 2: Constrained static relaxations of P-vacancy complex with P moved over various barriers. E™-^^ gives the height of the 
barrier for the P to move from a;*'' neighbour to the y*'', E^^ gives the binding energy of P at the a;"* neighbour site. The final 
panel shows the barrier for vacancy movement in pure iron. The barrier against motion of the vacancy-phosphorus complex is 
similar to the binding energy, and much lower than the barrier in pure Fe. Consequently, the P traps the vacancy, reducing 
overall vacancy diffusion. 
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FIG. 3: Mean Squared Displacement (MSD) against time for migration in a unit cell containing 2000 Fe and IP atom, initially 
located in a mixed dumbbell. Various temperatures were examined: at 600K and 920K the P remains attached to the interstitial 
throughout the period of the run. For 1200K, we observed rapid ID diffusion of the phosphorus atom via mixed interstitial (Fe 
does not diffuse in this phase), followed by a dissociation event and standard interstitial diffusion in pure Fe. The simulation 
box is large enough that the interstitial is not recaptured, and subsequent diffusion is for the iron atoms. 
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Potential 


Value (eV) 


Cutoffs (A) 


VFeFe{r) 


9734.2365892908C(r, 26, 26) /r 


0.0-1.0 




+ exp(7.4122709384068 - 0.64180690713367r - 2.6043547961722r^ + 0.6262539393123r^) 


1.0-2.05 




-27.444805994228(2.2 - rf 


2.05-2.2 




+15.738054058489(2.3 - rf 


2.05-2.3 




+2.2077118733936(2.4 - rf 


2.05-2.4 




-2.4989799053251(2.5 -r)3 


2.05-2.5 




+4.2099676494795(2.6 - rf 


2.05-2.6 




-0.77361294129713(2.7 -0=* 


2.05-2.7 




+0.80656414937789(2.8 - rf 


2.05-2.8 




-2.3194358924605(3.0 - rf 


2.05-3.0 




+2.6577406128280(3.3 - rf 


2.05-3.3 




-1.0260416933564(3.7 -0=* 


2.05.3.7 




+0.35018615891957(4.2 - rf 


2.05-4.2 




-0.058531821042271(4.7 - rf 


2.05-4.7 




-0.0030458824556234(5.3 - rf 


2.05-5.3 


VFep{r) 


(5615.9057245908/r)^(r, 26, 15) 


0.0-1.0 




+exp(10.76185442488 - 10.004045788895r + 4.9854254472397r2 - 1.2599788569372r3) 


1.0-2.0 




-3.3136605743629(5.3 - rf 


2.0-5.3 




+ 12.62523819360(5.3 - rf 


2.0-5.3 




-20.361693308072(5.3 - rf 


2.0-5.3 




+17.629292543942(5.3 - rf 


2.0-5.3 




-8.8120728047659(5.3 - rf 


2.0-5.3 




+2.5494288609989(5.3 - rf 


2.0-5.3 




-0.39698390783403(5.3 - rf° 


2.0-5.3 




+0.025779015833433(5.3 - r)" 


2.0-5.3 


Vpp{r) 


(3239.9456103409/r)C(r, 15, 15) 


0.0-0.9 




+e2;p(9.9382842499617 - 8.5637164272526r + 3.451962728599r2 - 0.61453831350215r^) 


0.9-2.5 




-0.078293794709143(5.3 - rf 


2.5-5.3 




+0.037557214911646(5.3 - rf 


2.5-5.3 


4lFeFe{r) 


11.686859407970(2.4 - rf 


0.0-2.4 




-0.01471074009883(3.2 - rf 


0.0-3.2 




+0.47193527075943(4.2 - rf 


0.0-4.2 


<l)Fep{r) 


<l,FeFe{21/2Af 


0.0-4.2 


<j)pp{r) 


<j}FeFei2l/2if 


0.0-4.2 


Fpeip) 


-yp- 6.7314115586063 x lO^^'p^ + 7.6514905604792 x W'^p^ 




Fp{p) 


-yp + 0.0011950274540243p^ 





TABLE V: Analytic form of the potentials. ^ is the short-range screening function[33 described in the text. As further work 
is done, we intend to continue the self-consistent fitting process to improve the potential. Users are encouraged to contact the 
authors with results and an evolving "best" set of parameters will be maintained online at http://homepages.ed.ac.uk/graeme 



